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NONLINEAR EVOLUTION OF HYDRODYNAMICAL SHEAR FLOWS IN TWO DIMENSIONS 
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ABSTRACT 

We examine how perturbed shear flows evolve in two-dimensional, incompressible, inviscid hydrody- 
namical fluids, with the ultimate goal of understanding the dynamics of accretion disks, as well as other 
astrophysical shear flows. To linear order, vorticity waves are swung around by the background shear, 
and their velocities are amplified transiently before decaying. It has been speculated that sufficiently 
amplified modes might couple nonlinearly, leading to turbulence. Here we show how nonlinear coupling 
occurs in two dimensions. This coupling is remarkably simple because it only lasts for a short time 
interval, when one of the coupled modes is in mid-swing, i.e., when its phascfronts are aligned with the 
radial direction. We focus on the interaction between an axisymmetric mode and a swinging mode. We 
find that all axisymmetric modes, regardless of how small in amplitude, are unstable when they interact 
with swinging modes that have sufficiently large azimuthal wavelength. Quantitatively, the criterion 
for instability is that \ky ^^vi / kx .aid\ ^ I'^/'Zl? i-^-, that the ratio of wavenumbers (swinging azimuthal 
wavenumber to axisymmetric radial wavenumber) is less than the ratio of the perturbed vorticity to 
the background vorticity. If this is the case, then when the swinging mode is in mid-swing it couples 
with the axisymmetric mode to produce a new leading swinging mode that has larger vorticity than 
itself; this new mode in turn produces an even larger leading mode, etc. We explain how this shear 
(or Kelvin-Hclmholtz) instability operates in real-space as well. The instability occurs whenever the 
momentum transported by an energy conserving perturbation has the sign required for it to diminish 
the background shear; only when this occurs can energy be extracted from the mean flow and hence 
added to the perturbation. For an accretion disk, this means that the instability transports angular 
momentum outwards while it operates. We verify all our conclusions in detail with full hydrodynamical 
simulations, done with a pseudospectral method in a shearing box. Simulations of the instability form 
vortices whose boundaries become highly convoluted. Whether this nonlinear instability plays a role in 
accretion disks is an interesting possibility. 

Subject headings: accretion, accretion disks — instabilities — solar system: formation — turbulence 



1. INTRODUCTION 

Matter accretes onto a wide variety of objects, such as 
young stars, black holes, and white dwarfs, through ac- 
cretion disks. For matter in an accretion disk to fall in- 
wards, angular momentum must be transported outwards. 
In many accretion disks, it is turbulence and the resulting 
turbulent viscosity that is responsible for angular momen- 
tum transport. Therefore, if one wishes to understand 
accretion disks, one must first address how disks become 
turbulent, and how much turbulence they sustain. In 
highly ionized disks, magnetic fields can trigger turbulence 
via t he magnetorotational instability (jBalbus fc HawlevI 
Il998() . However in neutral (hydrodynamical) disks, the 
situation is unclear. It remains an open question whether 
neutral Keplerian disks are turbulent, or whether they can 
remain laminar despite enormous Reynolds numbers. Un- 
til it is answered, the accretion of nea rly neutral disks , 
such as those aroun d young stars fe.g.. JSano et ani200Cl[ ) 
or dwarf novae fe.g..lGammie fc Menoulll99 8V will remain 
myst eries. Simulations (Hawle v et al, [199 9: S hen et al.l 
|2006[) and experiments (ji et al. 1120061) both point to the 
answer that incompressible hydro disks are laminar — or 
at least that any turbulence in them would not be strong 
enough to act as the agent of angular momentum transport 
(jLesur fc Longaretti 2005) . But the evidence is not con- 
clusive because the Reynolds numbers in disks are much 



larger than those accessible to computers or experiments. 

Even if one takes the view that numerical and exper- 
imental evidence rules out self-sustaining incompressible 
turbulence in Keplerian disks, it remains important to un- 
derstand how such shear flows work. In our view, this is 
an essential first step before more complicated effects — 
such as baroclinic effects, convection, stirring by planets, 
stirring by dust, or stirring by a magnetically active layer 
above the midplane — can be understood. It could be that 
one of these effects is responsible for turbulence in disks. 
In the present paper, we describe in detail the dynam- 
ics of perturbations in incompressible shear flows in two 
dimensions (in the plane of the disk). Even though 2D 
flows do not lead to turbulence, they exhibit a remark- 
ably rich dynamics. Of course, incompressible 2D shear 
flows have be en the sub ject of an enormous amount of re- 
search (e.g., iDrazin fc Rcid 2004). But our methods are 
different: we calculate the nonlinear coupling between lin- 
ear waves. One of the beneflts of this approach is that it 
is straightforward to generalize it to three dimensions, as 
well as to include effects such compressibility and magne- 
tohydrodynamics. 

Recently, the evolution of linear waves has been the 
focus of much attention in the astrophysical literature. 
Since the velocity field of a linear swinging wave can 
be amplified by very large factors during its swing, 
it has been proposed that sufficiently amplified modes 
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might cou ple nonlinearly , leading ultimately to turbu- 
lence (e.g. 



loannou fc Kakouris 2001; Chagelishvili et al 



2003HUmurhan fc Regevll2004l: lYeckdl2004i; lAfshordi et al 



20051) . However, how this coupHng might occur has not 



been discussed (though see lMukhopadhvavll2006D . nor has 
the amplitude needed to trigger nonlinear effects been 
quantifie d. The present paper beg ins to rectify these short- 
comings. iBalbus fc HawlevI ((2006') show that linear swing- 
ing waves are exact solutions of the nonlinear equations of 
motion. They contend that this argues against transient 
amplification as a route to turbulence in unmagnetized 
disks. However, we show in the present paper that non- 
linear coupling between two different waves can lead to 
interesting dynamics even when the individual waves are 
exact nonlinear solutions. 

Numerical simulations of 2D shear fiows starting from 
random in itial conditions settle into a distinctive banded 
struc t ure (lUmurhan fc RegevI |2004 iJohnson fc Gammi^ 



120051 : IShen et al 20061 ) . Roughly speaking, bands where 



to > are interspersed with bands where w < (w is the 
perturbed vorticity, eq. [7] ) . Bands where lu has the same 
sign as the background vorticity contain a single vortex; 
in bands with the opposite sign, uj is smooth and there are 
no vortices. One of the goals of this paper is to explain 
why this is a natural outcome of random initial conditions. 

1.1. Organization 

We introduce the equations of motion in fj2l In §fj3]- 
[HI the heart of this paper, we describe and simulate the 
nonlinear coupling between modes, focusing on the shear 
instability that results from the coupling between swing- 
ing waves and axisymmetric ones. In SjH which can be 
read independently of §[210 we describe the instability in 
real space. We give both a dynamical explanation and one 
based on momentum and energy arguments. We also sim- 
ulate the fully nonlinear outcome of the instability. We 
conclude in ^JT) In the Appendix, we describe the pseu- 
dospectral code that we use to run simulations in §^5]l5l 

2. EQUATIONS OF MOTION 

An unperturbed Keplerian disk has angular velocity pro- 
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We write the fluid equations in a refer- 



ence frame rotating at constant angular speed fig = ^^('^0)7 
where rp is a fiducial radius, and replace the radial and az- 
imuthal coordinates r, 9 with Cartesian coordinates, x = 
r ~rQ, y = ro9. On lengthscales much smaller than tq, the 
"shearing box" equations of motion for an incompressible 
fluid are 

dtv + v ■Vv = -2noZXv + 2qnoxx-V{P/p) , (1) 
V-v = 0, (2) 
where v is the velocity field in the rotating frame and 
dfl 



q 



d\nr 



(3) 



We use standard Cartesian vectorial notation, e.g., v = 
v,j.x + Vyif + VzZ, etc. The first term on the right-hand 
side of equation ^ is the Coriolis force, and the second 
is what remains after adding the centrifugal and gravita- 
tional forces. 



Decomposing the fluid velocity into 



-qxy 



(4) 



where the first term is the shear flow of the unperturbed 
disk, yields 



[dt - qxdy) u + ■ 







(2ilo - q)u^ 



y - y{pm 

(6) 



We shall solve these equations in two dimensions, in the 
X — y plane. It is simpler, though, to work with the curl of 
equation ([5]) , which may be expressed in terms of vorticity 
of u, 



as follows: 



Lj = z • (V X u) = dxUy — dyi 



{dt — qxdy)Lo + u • VtJ = 



(7) 



(8) 



with u given by the inverse of equation ([7]) : 



z X VV" 



(9) 

This paper is de- 



Equations (I8])-([9|) form a complete set. 
voted to solving them. 

Equation ^ shows that to is advected by the total ve- 
locity field —qxy + u. Vorticity is locally conserved in 
two dimensions (in the absence of dissipation); to is nei- 
ther created nor destroyed, but is simply shuffled around 
by the velocity field field that it induces via equation 
Therefore in 2D investigations one is forced to specify an 
initial vorticity field, and then to see how it is shuffles itself 
around. In three dimensions, vorticity might be created 
by turbulence, or by stirring by planets, or by convection. 
But the creation of vorticity is beyond the scope of this 
paper. 

The vorticity of the unperturbed flow in the rotating 
fram^l is 



z ■ Vxi-qxy) = -q ; 



(10) 



q sets a scale for the vorticity against which the strength 
of the perturbation lo may be measured. Throughout this 
paper, we restrict the perturbed vorticity to be less than 
the background vorticitjd, 



< 



(11) 



It is conceivable that hydrodynamic disks are so turbulent 
that they violate this inequality. But if they reach such a 
state via a nonlinear instability, then the instability would 
have had to act while the inequality held. Because of the 
above inequality, it is tempting to drop the nonlinear term 
u • Vw from equation ([5]). However, we shall show that 
even when ^ q the nonlinear term is important pro- 
vided that u; is sufficiently elongated in the streamwise 
direction. While we may indeed neglect nonlinear advec- 
tion in the y-direction (uydyUj), the term UxdxOJ advects 
fluid in a direction orthogonal to the background shear. 
The net displacement in x can be significant if Ux acts 
coherently over a large range in y. 



In the non-rotating frame, the unperturbed disk has tot al v elocity field utot = rflO, and the corresponding vorticity is z ■ VXDtot = 
dQ/dlnr + 2f2 = —q + 2Qo at rg. This differs from equation ||10|I (by 2^0) because of the change in reference frame. 

^ This restriction precludes us from seeing the shear (or Kelvin-Helmholtz) instability discovered bv IShen et al.l l|2006l) . which requires \uj\ > g. 
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2.1. Equivalence with Non-rotating Linear Shear Flows 

Equations dS])-® also describe the dynamics of non- 
rotating incompressible 2D linear shear flows. To see this, 
one may discard the first two terms on the right-hand side 
of equation ([1]) and make the decomposition of equation 
([4]), where the unperturbed shear rate q may now be set 
to any arbitrary constant. Equation ([5]) is then repro- 
duced without the coriolis terms (= —2^}oZ X u). But 
the curl of this equation still yields equation ([5]), because 
in two dimensions the coriolis force is the gradient of a 
scalar, —2z X u = VV~^2ti;. Therefore the coriolis force 
can be absorbed into a redefinition of the pressure, and 
does not affect the dynamics. (The sum of centrifugal and 
gravitational forces 2qrtoX also does not contribute to the 
dynamics of u — it merely cancels the coriolis force due to 
the background shear —2i}oZX{—qxy). But this is true 
in 3D as well. In fact, this balance sets the unperturbed 
velocity field.) 

Rotating and non-rotating incompressible shear fiows 
are equivalent only in two dimensions. To illustrate how 
this equivalence breaks down in three dimensions, we ex- 
amine linear axisymmetric waves. In non-rotating shear 
flows (eq. [S] without the coriolis force) , such "waves" have 
dispersion relation uj — 0, and eigenfunction u oc y. They 
are static disturbances of the background flow that merely 
alter the shear flow's velocity profile. Rotating shear flows 
have dispersion relation = 2f7o(217o — g)fc^/fc^. There- 
fore two-dimensional modes (with kz = 0) have uj = 
whether or not the flow is rotating. But axisymmetric 
modes in rotating flows with fcz ^ have a non-zero fre- 
quency. When three dimensional motions are allowed, the 
coriolis force induces epicyclic-like oscillations. 

3. NONLINEAR MODE COUPLING EQUATION 

We solve equations ©-([Q]) in a box of size x Ly sub- 
ject to boundary conditions that are periodic in y, and 
shifted periodic in x (eq. |A5| ). Lagrangian coordinates 
{X,T) shear with the background flow: 

Y = {y + qtx) mod Ly (12) 

{X,T) = {x,t) (13) 

^dr^dt-qxdy . (14) 
The Fourier transform of uj with respect to Lagrangian 
coordinates is 



where 



Jo Lx Ly 



^^^^X)e-^K.X 

T T ^ ) 



211 Jq- 



for integer Jx, Jy, with inverse 



E 



iK ■ X 



(15) 
(16) 

(17) 



Capital letters {X,K,Jx,Jy) denote Lagrangian coordi- 
nates, and lower-case letters denote Eulerian ones. The 
transform of equation ([S]) gives the nonlinear mode cou- 
pling equation, 

di^J.,J. V- z- [k{K - K',T)xk{K 



where K' is related to (J^, J'y) via the analog of equation 
mi), and 

k{K,T) = K + qTKyX (19) 

is the wavevector with respect to x. Equivalently, since a 
single mode has the spatial dependence 



exp(ilf • X) = ex.p{ik{K,T)-x) , 



(20) 



the Fourier transform of V = {dx,dy) is obtained by 
making the substitution V — > ik(K,T). Hence equation 
(fT8|) can be read directly from equation ([8|) after replacing 
the Fourier transform of u with {—i/k'^)z X k times the 
Fourier transform of oj (eq. [9]). 

4. LINEAR EVOLUTION 

The linearized equation of motion is duj/dT — 0, which 
has the general solution is 

L0 = f{X,Y)^ f{x,v + qtxmoALy) (21) 

for any function /. A single vorticity wave can be written 

as 

u) = UJ cos {K • X) = UJ cos {k{K ,T)-x) , (22) 

where uj and K are constants. The velocity fluctuation 
due to this wave is 



u = z X VV UJ 
_ z X k{K,T) 
\HK,T)\^ ' 



'Sin (K -X) 



{^Ky,Kx + qTKy) 



wsin (K ■ X), 



(23) 
(24) 

(25) 



{Kx + qTKyY + Kl 

the amplitude of which depends explicitly on T . 

Equation ([22|) is not only a solution of the linearized 
equation of motion, but it is an exact solution of the non- 
linear equation as well: the nonlinear term u • Vo; van- 
ishes for a single wave because of incompressibility. 

4.1. Axisymmetric, Swinging, and Radial Waves 

Axisymmetric waves have ky — Ky — 0: their phase- 
fronts are aligned with the y-axis, which corresponds to 
the azimuthal direction in a Keplerian disk. See Figure 
[TJ Axisymmetric waves are independent of time not only 
in shearing coordinates, but in unsheared ones as well. 
"Swinging waves" have Ky 7^ 0. At early times (T 
—00), they are leading waves: the ratio kx/ky — > —00, so 
their phasefronts, though nearly axisymmetric, are tilted 
slightly into the first and third quadrants of the x — y 
plane. As time evolves, the phasefronts swing first through 
the radial direction — parallel to the a;-axis — and then, as 
T — > -|-oo they approach re-alignment with the azimuthal 
direction as a highly trailing wave. They are exactly radial 
at time 

r,ad(K) ^ -4^ , (26) 



qKy ' 



dT 



terms of which 



|fc(K',T)|2 
z- [K X K'] 



kx = KyqiT - r,ad) . (27) 

We refer to a swinging wave at time Tjad as a "radial" 
(Ig^ave, since its phasefronts are then aligned with the ra- 
dial direction. 
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Simulation Parame ters 

n =30 
n;=10 

Lj=87il00 

timestep= 1/150 
S=0.6i 
J.= 10 
I,= -l 

K,sZ7iJ,/L^=10 
■• 2;iJ/L,=-0.01 




Fig. 1. — Linear Evolution in fc-space: Only the half of fc- 
space with fe^: > is plotted, the other half being redundant because 
a mode with wavevector — fc is the complex conjugate of a mode with 
wavevector fc. Two modes are shown. The axisymmetric mode is the 
white square at {kj:,ky) = (1,0); it does not move in fc-space. The 
nearby rectangle depicts its contours of constant vorticity in x — y 
space, and shows phasefronts aligned with the background shear, 
parallel to y. The second mode is a swinging mode. It evolves in 
time from the lower-right corner to the upper-right one. Rectangles 
again depict phasefronts at the corresponding positions in fe-space. 
The dashed arrow at kx = can be thought of as an instantaneous 
jump in ky-. as a mode with {kx, ky) = (0, —1) continues off the plot 
to negative k^, its complex conjugate appears at (0, 1). 



4.2. Numerical Simulation: a Single Swinging Wave 

Figure [2] shows the evolution of and Uy of a swinging 
wave. The points show the output from the pseudospectral 
code described in the Appendix, and the curves through 
the points are given by equation ([^5]) . See the Figure's 
caption for details. Figure [3] shows contours of constant 
vorticity at three times in this simulation. 



5. NONLINEAR EVOLUTION 

The nonlinear coupling coefficient between two modes is 
time-dependent (eq. [H]), which at first sight makes the 
mode coupling problem look forbidding. But the situation 
is not quite so dire. The coupling coefficient depends on 
time only through its denominator 



\k{K',T)\-'=K''(l 



f{T- 



(28) 



which is smallest at the time when either of the modes has 
radial phasefronts {k^ = 0). At much earlier or later times 
the coupling oc . Therefore integrated over time, most 
of the coupling occurs around the time that = 0. Be- 
cause mode couplings occur at essentially discrete times, 
the analysis is greatly simplified, as we presently slow. 



Fig. 2. — Evolution of a Single Swinging Wave: At T = 0, 

a vorticity wave (eq. 1221 ) is initialized with parameters as shown 
in the figure inset. Jx is the value of = KxLx/2it of the only 
non-zero mode; and similarly for Jy; n^, riy are the number of grid- 
points in the simulation, and time is measured in units of i^g^ ■ 
Points show output from pseudospectral code, properly normalized, 
i.e., they show 2Im(iiio,_i), where uj^^j^ is the Fourier transform 
(as in eq. 1151 ). Lines through the points show the amplitude of 
equation I I25I I: ibKy/k"^ in the top panel and Qk^/k'^ in the bot- 
tom. At time Tj-ad = —Kx/qKy = 667, phasefronts cross through 
the radial direction. In the pseudospectral code, dealiasing with the 
2/3 rule sets itio,_i to zero when this mode has \kx\ > fia;/3 = 10 
feg. IA2l1 ). This occurs at time 2f)/q\Ky\ = 1333. At a later time, 
when \kx\ = n^/'i = 15 (which occurs at time 2b/q\Ky\ = 1667) the 
kx of this mode is mapped from fcj; = — 15 to kx = 15 via the mod- 
ulus function in equation IIA17I I. So this mode once again becomes 
a leading mode. However, since the amplitude of this mode was set 
to zero by dealiasing, it remains equal to zero. 



5.1. Swinging Waves and a Single Axisymmetric Wave 

To begin our investigation of the nonlinear coupling be- 
tween waves, we consider the evolution of small swinging 
waves in the presence of a much larger axisymmetric wave. 
We set 



UO = 2cJi,oCOs(X) + UJswing,{X ,T) , 



(29) 



with l^swingl ^ ^ li and oji q a real- valued con- 

stant; the factor of 2 multiplying wi.o has been inserted 
for consistency with equation (fT5|) . We have chosen the 
length unit so that the axisymmetric mode has — 1. 
Substituting equation ((29)) into equation ([8]) gives 



dT 



-2tJi,osin(X) (V^^ -f 1) dyuj, 



■y^swing 1 



(30) 



after dropping itswing-Vwswmg- 
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Fig. 3. — Evolution of a Single Swinging Wave: Contours 
of constant vorticity at three times from the simulation described in 
Figure [2] showing, respectively, leading, radial, and trailing phase- 
fronts. 

Equation (j30p is linear in Wswing, with a spatially vari- 
able cocfBcicnt, sin(X). A Fourier mode of Wswing with 
wavevector {Kx,Ky) couples with sin(X) to produce two 
modes with wavevectors {K^ ± 1, Ky). Hence a solution of 
equation ([50]) i^ 



E 

— — oo 



2ujj^ (T) cos {J^X 



KyY) 



(31) 



Since we take Ky of the swinging modes as fixed (with 
Ky = —Ky < 0), we suppress the corresponding subscript 
Jy from wj^, ; mode amplitudes with a single subscript all 
label modes with the same Jy ^ 0, where Ky = 2'Kjy/Ly. 
Substituting the above Wgwing into equation (|5D|) . 



dujj 



dT 



= -K 



""Jx-l-l 



where 



and 



qTKy 



q 



rad , J_ 



qKy 



(32) 

(33) 
(34) 

(35) 



(eq. [26]). In writing equation ([32|) . we kept only the V ^ 
term from the right-hand side of equation ([50]) . To include 
the full -I- 1, one should replace the kj^^-^ in equation 
(|32p with k^\i — 1; we show below that the omitted terms 
are small (eq. [46]). The term is due to the advec- 
tion of the axisymmetric wave by the swinging ones, i.e., it 



is — Mswing" V2wi,o cos (X). The dropped term represents 
advection of the swinging waves by the axisymmetric one. 

If Ky <^ 1 , then the squared wavevectors that appear in 
the denominator in equation (j32p reach very small values 
when the corresponding mode is radial. At this time, the 
mode strongly influences its two neighboring modes. We 
shall see that we can often approximate the interaction 
between different swinging waves as occurring solely when 
one of the waves passes through the radial direction. 



As an example, at T 
single leading wave: 



we initialize uig 



with a 



^swing 



(X, T = 0) = 2e cos{X - KyY) 



i.e., with oji — e, and i^j^_= for ^ 1. We choose 



0<Ky^l 



and follow the evolution governed by equation 
Figure [H 



(36) 

(37) 

(38) 
See 



0-- 



^0 = Oi^i)= -3i 

□ 



-if. 



■W2= ex, 



Fig. 4. — Early-Time Evolution in fc-space: The top panel 
shows four modes between the times and 3/AqKy. For 
example, the square at {kx,ky) = (l,~Ky) represents uji at 
time T = 0; the connecting arrow shows this mode's tra- 
jectory in fc-space until time Z/AqKy. The middle panel 
continues the evolution through time Tj-^^ i = 1/qKy, when 
Lui is radial, and cjo and L02 interact strongly with the radial 
and axisymmetric modes. At this time, 0^2 changes from 
to €Xi, and loq changes from to —exi- The bottom panel 
follows the evolution through time Tj^d 2 = '^/qKv 



This is clearly not the complete solution. But each "chain" of modes represented by equation (I31I I. with fixed Ky and 



< Jj; 



< 00, 

evolves independently of others. For the complete solution, one should also sum over different Ky, and include chains whose kx{T = 0) are 
not integer multiples of the axisymmetric mode's K^', these chains are also described by equation 11311 . except that their time origin (when 
Kx = kx) is offset from T = 0. Note also that for simplicity we have not included arbitrary phases in the arguments of the cosines (i.e., we 
take mode amplitudes to be real- valued) . It is trivial to include phases; we do so in t|5.4l 
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At first, the only non-zero swinging mode is tui. This 
mode influences only two others: lj2 and cuq. So at early 
times, 0J2 approximately satisfies 
duj2 



Ky 1 + 



(39) 
(40) 



Integrating this from time T = through T,., 
l/{qKy), we see that most of the integral comes from 
within a time ~ 1/q oi Trad.i, when the phasefronts of 
bJi are within around ±45° to the radial direction. So 
we can extend the integration to the domain —00 < 
T < 00 to obtain shortly after time Trad.i, i-e. at time 

(1 + qKy)Trad,l, 

W2 = ujixi = exi (41) 

where 

(42) 



Xi 



qK, 



is the dimensionless parameter that controls the linear evo- 
lution of Ci^swing- The behavior of luq at time Trad.i is sim- 
ilar to 0^2, with = ^eXi shortly after Trad,i- 

In the top panel of Figure 21 the modes initially evolve 
linearly, with constant amplitude and = — qTKy. 
Before time Tiad,!: both wi and lj2 are leading modes 
(their k^/ky < 0), and is trailing. But when wi be- 
comes radial and its k^ = 0, the amplitudes of modes 
loq and UJ2 change abruptly (Figure [H middle panel). 
Qualitatively, the velocity field of the uji mode when 
it becomes radial is u'^'^'^ — —x{2e/Ky)sm{Kyy) (eq. 
P5]l. Since it takes a time ~ 1/g to swing through 
the radial direction, the corresponding displacement field 
is ^^'^'^ ~ u^^'^/q, which advects the axisymmetric 
mode, changing it by an amount — ^''''^•V2cji_o cos(X) ~ 
—2e{uJifi/qKy) (cos(a; — Kyu) — cos{x + Kyu)). Hence 
the UJ2 mode, which has the spatial dependence 
cos {x — Kyyj at this time, changes its amplitude by ~ 
—e{ijJifi/qKy), as in equation ([1T|) : similarly, ujq is changed 
by an equal and opposite amount. 

After time Tiad.i, the next time of interest is Trad, 2 = 
IjqKy^ when the 102 mode is radial with 102 — eXi 
(Figure SI bottom panel). The evolution of ^3 around 
this time is given by the analog of equation (|39p. i.e., 
dhj'i/dT — —KyUJi,QUj2/k2- Integrating through time Tiad,2 
yields 1^3 = UJ2X1 = ^Xi shortly after Ti.ad,2- Similarly, ui 
changes by an equal and opposite amount. 

Extrapolating to later times, it is clear that the ampli- 
tude of the mode that crosses through the radial direction 
is 



at time T^. 



J X 

qKy 



(43) 



As time progresses, the amplitudes of the radial waves 
grow exponentially if |xi| > 1; conversely, they decay ex- 
ponentially if Ixil < 1- At marginal stability, 

= ^ . (44) 

marginal stability Q 



: 1 for the axisymmetric mode then 
Since ~ n/Ky, equation (|43l) 



(If we do not set 

K = TrlXjj^axit^l.o/?!)- 

shows 

(X exp (|wi,o|T'rad,J,7r(^y/K)ln(K/^y)) (45) 
and hence the growth rate is very small both when Ky k 
and when Ky 0; it is fastest for Ky = n/e = 0.368k. 



5.2. Numerical Simulation 

Figures [5] and [6] show results from two pseudospectral 
simulations. At T = 0, the vorticity field was given by 
equations ^ and i.e., uj{T = 0) = 2wi,oCos(X) + 

2ecos — KyYy One of the simulations had Xi = ^1-1, 
and the other had xi = —0.9. In Figure El we plot the 
amplitudes of the first eight modes that pass through the 
radial direction in the xi ~ —1-1 simulation. Most of 
the time the amplitudes remain constant. Only at times 
?iad,j^ do the amplitudes of modes Ja;±l change abruptly; 
see also Figure HI 
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Fig. 5. — Swinging Waves in the Presence of a Large Ax- 
isymmetric Wave: Pseudospectral simulation of the fully 
nonlinear equations of motion. The vorticity field was ini- 
tially uj{T = 0) = 2a;i,o cos(X) -f 2e cos (X - i^'yy) (eqs. [251 
and I36h : parameters are listed in inset. Curves show 
LUj^ = Re(a;j^ Dealiasing is described in flAl 

The points in Figure [S] show the amplitudes of radial 
modes at each time that a mode is precisely radial. The 
line through the points is the theoretical prediction. (Of 
course, only at the discrete times T^s^d.j^, for integer J^, is 
there a radial mode; the theoretical prediction is plotted 
as a line for clarity.) The theoretical prediction is given by 
equation (|43p. but with a slightly modified xi- In partic- 
ular, in deriving equation (I43|) . we made two approxima- 
tions that are applicable as Ky — > 0; we shall now include 
the 0{Ky) corrections. First, we integrated equation ((39)) 
from —00 to 00, whereas we should have integrated from 
~ to ~ 2/ qKy] second, instead of uJi/kf in equation ([5^ . 
we should have wi(l/A:f — 1). Both corrections can be ab- 
sorbed into a redefinition of the dimensionless parameter. 



Xl = Xl • (1 Ky) 

TT 



(46) 



Xi should be inserted in equation (j43p in place of xi- As 
long as Ky ^ 1, this correction is small. Nonetheless, 
accumulated corrections can be substantial if xi is expo- 
nentiated many times. For example, in the top panel of 
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Figure [HI using xi = —1-1 instead of xi = —1.086 would 
overpredict the amplitude of the last plotted point by a 
factor of 2. 




Fig. 6. — Amplitudes of Radial Modes: Points in the top 
panel show results from the pseudospectral simulation de- 
scribed in Figure [Sj extended to later times. Each point 
shows the amplitude of the mode that is radial at that 
time. The line through the points shows the theoretical 
prediction e|xi|'^""^~^, as in eg. 1431 . but with xi = —1.086 
in place of xi = —1-1; see eq. | 46 |. Bottom panel shows 
results from a simulation identical to top panel's, except 
with a;i,o = 0.0043 => xi = -0.9 and e = lO^". Note that 
the last point in both panels has Jx = 58, even though the 
simulations only have nx/2 = 15 modes in the x-direction: 
modes are "recycled" by the r emap ping that is performed 
by the modulus operator (eq. IA17n . 

5.3. Quasi- Eulerian Notation 

Analysis of the coupling between Fourier modes is sim- 
pler when a quasi-Eulerian notation is used; we introduce 
this notation here. Instead of uj^j^, where the subscripts 
label the mode's Lagrangian wavevector K, we label a 
mode's amplitude 

^j^-Jy (47) 

when the mode's Eulerian wavevector k{K,T) satisfies 



(48) 



J" / 
^ / 

-1 



,0.1 _/. ,0.-1 I I 

i^->0 o — 

, ,i.O 



Q Q 



3 



Fig. 7. — Quasi-Eulerian Notation: Squares represent 
the wavevectors fc of a number of modes at two successive 
times, T = T^ad - l/iqRy and T = T^ad + l/iqKy, where T^ad 
is the radial time of the mode labelled u)^'^^ in the Figure. 
Modes are labelled by superscripts jx,jy representing the 
nearest kx,ky. Axisymmetric modes have ky = 0. Because 
their fc does not evolve, their quasi-Eulerian indices are the 
same as their Lagrangian ones. 

5.4. Many Axisymmetric Waves 

We have shown above that a single axisymmetric wave 
is unstable to interactions with swinging waves whose Ky 
are sufficiently small that \xi\ > 1- In the present sub- 
section, we show how to generalize this result to consider 
many axisymmetric waves, and so re-derive standard re- 
sults regarding the instability of shear flows with inflection 
points in the limit that ^ g and that the vorticity is 
elongated in the streamwise direction. 

We generalize equation (^5]) to 

= Waxi(a;) -I- Wswing(^, T) , (49) 

where ojaxi is an arbitrary axisymmetric function, 



(50) 



where and jy are any integers. Since k^. is time depen- 
dent, equation (|48p only applies at discrete times. But we 
shall also label a mode's amplitude with uj^^^iv at times 
shortly before, and shortly after, its fc satisfies equation 
((48|); see Figure [71 (To be more precise, one might say 
that a mode whose increases in time has its label 
switched from oj-'^ to uj^'^+^'iy when its k^ passes through 
{2tt/ Lx){jx + 1/2); but this increased precision is unnec- 
essary.) 

^ We are assuming that the amplitude of a swinging mode hardly chang 
This assumption follows from the requirement that Ky be much smaller 



with u!^'"'^ arbitrary constants and = 2tt. We seek so- 
lutions to linear order in Wswing- No new values of ky can 
be generated by mode couplings because swinging modes 
must couple with axisymmetric modes, which have ky = 0. 
We take Wswing to have ky — iKy, where Ky > 0, and Ky 
is much smaller than the typical k^ of Waxi- Equation p8p 
gives the equations of motion in Lagrangian notation. We 
switch to quasi-Eulerian notation by replacing subscripts 
with superscripts, and capitalized indices with lower-case 
ones. As in for small Ky most of the coupling be- 

tween modes occurs when one of the modes is radial, i.e., 
when the term \k{K' , T)p in the denominator of equation 
(|18|) nearly vanishes. So we take ujji^^j^ in this equation to 
be the nearly radial mode w"'^^ (setting Ly — 2'K/Ky so 
that jy = ±1), and, from equations (fT8)) and (p7| . 

= -c."-'ic.-."i^ -i , (51) 

dT Kyl + q^{T-T,^if ' ^ ' 

where Trad is the time when lo'^'^^ is radial. Equation 



uj) generalizes the first term in equation (|32p to arbi- 
trary jx and complex mode amplitudes. We integrate from 

:es as it crosses through the radial direction, i.e., for —Ky ^kx ^ Ky. 
than the typical kx of ojaxi . 
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T = — oo to +00 because most of the integral comes from 
within a time ~ 1 of Tj-ad The change in uj^'= ■ ^ ^ at this 
time is 



where 



0,- 



qK, 



(52) 
(53) 



y 

0,-1 



The mode that is labelled w at time Trad was, at the 
earlier time T^ad — ^/iKy, labelled uj^'^^. At that ear- 
lier time, its change due to the then-radial mode was 



Xi^ ' I 
we conclude 



-l/qKy 



Extrapolating to still earlier times. 



T=J/qKy 



E 



XjL 



(54) 



where J is an integer. Solutions of this difference equation 
have the form 



.0, 



— constant x z"^ , (55) 

T=J/qKy 

with z determined by the roots of the characteristic equa- 
tion 



1 = 



00 



(56) 



This is the dispersion relation; unstable modes have \z\ > 
1. The eigenfunction is, generalizing equation ([54|) FI 



T=J/qKy 



00 



constant x z"" 



T={J+j^ 



E 



-(57) 

-r.)/QKy 



At marginal stability. 



-i/3 



z = e 

for a real-valued /?, and the dispersion relation (eq. 
is 



1 = 



1 



2qKy 



dx 



dx 



x=l3 



(59) 



(60) 



- 13 2qKy 

after changing the sum over j'^ to an integral, and ex- 
tending the box-size L^; to 00. The imaginary part of the 
dispersion relation shows that /3 must be chosen at a zero 
of dhJayii/ dx] the real part gives the marginally stable Ky: 



Ky 
1 

Yq 



marginal stability 

duj^y,i/dx 



dx- 



/3 



/3— zero of dto-^y^il dx 



(61) 
(62) 



stability criterion. Equation ([5^ is derived bv lGilll (jl965[ ). 
although using an entirely different method. That shear 
flows can become unstable only when the velocity field has 
an inflection point (i.e., when dujpi-x\ /dx = 0) is known a s 
Rayleigh's inflection point theorem (|Drazin fc Reidl[2004[ ). 
For an unstable mode, 

a > , (63) 



z = e"-"^ 



and the dispersion relation is 
1 



1 = 



2qKy 



dx 



diosxi/dx 
t: — (3 — ia 



(64) 



which reduces to equation ([60]) as a ^ 0. 

Since the eigenfunction is a convolution in Fourier space 
(eq. [58]), it is simply expressed in real-space: 



^aqTKy ^^axi 

dx 



Re const- 



,-iKy(y+f3qTKy) 

l + i{x- P)/a 



(65) 

where the constant is an arbitrary complex number, and 
we have set J — qTKy. A more transparent way of writing 
the eigenfunction is to combine it with cjaxi 

LO = Waxi + t^swing ^axi(a; + ^{x, V, T)) , (66) 

where the displacement field ^ is given by 



C(x,y,T) =e(0, 2/0,0) X 



Re 



1 



1 + ix/a 



'iKy(y-yo) 



I/O is an arbitrary constant, and we have set /3 = 0, which 
can be done by changing the origin of the x-axis. 

The exponential growth rate is 7 = aqKy. Note that 
a has the physical interpretation that fluid at x = ±a is 
advected by the background flow ~qxy a distance 1/ Ky in 
the growth time I/7. Fluid at \x\ 3> a is advected many 
wavelengths in a growth time; this phase mixing implies 
that the eigenfunction is cut off at |x| ^ 

5.5. Numerical Simulation 

We use the pseudospectral code to simulate the equa- 
tions of motion with Waxi given by a Gaussian profile, 

Waxi(a;) = -wexp (-a;^/s^) , (68) 

with sKy <C 1 and Co > Q. This profile is somewhat unre- 
alistic because it does not integrate to zero. Hence it pro- 
duces a velocity field Uy that is infinite in extent. Nonethe- 
less, we show in i j6.5l ttiat profiles that do not integrate to 
zero behave similarly to those that do. 

The marginally stable wavenumber is (eq. [62]) 

K = V¥— . (69) 
qs 

The imaginary and real parts of the dispersion relation 
(eq. [53]) imply 



1 : 







qKyS^ 



dx- 



1 1 2 

'X j s 



In gSTTl we had cjaxi ~ 2ciJi_ocosx; inserting this into the 
above integral gives k = 7r|cL;i_o|/(?, as in equation (I44|) . For 
an arbitrary profile of Waxi, if the right-hand side of equa- 
tion ([5^ is negative for each /3 that is a zero of diOaxi/dx, 
then that profile is stable to all small perturbations. Con- 
versely, if the right-hand side is positive for any /?, then 
there are unstable modes with Ky < k. This is the general 

^ To be more precise, the left-hand side of equation II57I I should be evaluated not precisely at time J/qKy, but slightly sooner — before this 
swinging mode has been altered by the radial mode. But this distinction is unimportant when the typical of io^i is much larger than unity, 
in which case may be treated as a continuous variable. 

^ Although ojs^ving is cut off at |x'| > a, awing extends further, to \x\ ~ ±1/Ky; see JS] 



^ f 1 _ V^-e("/^)'erfc(a/s) 

Ky \ S 



(70) 
(71) 

(72) 
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Solving the latter equation for a gives the exponential 
growth rate 



K 

7 = aqKy = u}f{—^) 

K 



(73) 



where f{Ky/K) is plotted in the top panel of Figure[8l The 
three points in that panel show results from numerical sim- 
ulations; see caption for details. The fastest growing mode 
has Ky/n = 0.299 and growth rate = 0.435(1'. The bot- 
tom panel shows that when our assumption sKy <C 1 is no 
longer valid, the growth rate differs from equation ([75]) . 
In Fourier-space, the eigenfunction at a fixed time is 



constant x e 
-a 



1 - ^/^■ 



-^e(aA+sjW2)^erfc(a/5 
s 



sjW2))74) 



after changing the sum in equation (|58p to an integral. In 
Figure ini we show that this agrees with the output from a 
numerical simulation. 

Figure [TO] shows contours of constant oj at two times, 
and shows that theory agrees with simulation. At later 
times than those shown, the vorticity tends to wrap up 
into a vortex, with large spatial variations m. lo. To simu- 
late this with the pseudospectral code, we must introduce 
an explicit viscosity, which smooths out the variations in 
oj. But the highly nonlinear state can be simulated with 
no viscosity with a Lagrangian code; we shall do so in 
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Fig. 8. — Exponential Growth Rate: Curve in top panel 
shows growth rate obtained by solving equation (|72|l for a/s 
with Mathematica. Points show results from three simula- 
tions with varying Ky, but otherwise identical. In the sim- 
ulations, tiJaxi was initialized as ~Q (cxp{— x^/s^) — ,Jtts/ Lx) , 
which differs from equation (I68II by the addition of a 
constant term that eliminates the zero-frequency compo- 
nent, i^swing was initialized as ecos{Kyy), with e = 10^''. 
The growth rate was extracted from each simulation by 
plotting the amplitudes of the radial modes versus time 
(as in Figure [SJ, and fitting with an exponential. For 
the fastest growing mode, Ky/K = 0.299. In the bottom 
panel, we show results from six simulations with fixed 
Ry/n = 0.299 => s Ky/u ) = 0.353, and fixed s, but varying Ky 
and il). Equation (I73D predicts that the growth rate should 
be equal to 0.435 iD, independent of sKy. Numerical results 
(points in figure) agree with this prediction for sKy < 0.01. 
The theory leading to equation (I73D is only applicable for 

SKy < 1. 
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Fig. 9. — Vorticity in Fourier-Space: Points show output 
from the simulation described in FigurelSlthat had u> = 0.01, 
s = 0.1, Ky/ft = 0.299, and hence a = 0.082. The vortic- 
ity shown was output at time T = 1208. Top and bottom 
panels depict the eigenfunction. By the reality condition 
(^jx,i = (w^J^'^^)*. The curve through the points in these 
two panels is equation (I74D . with the normalization chosen 
to fit the points. The agreement between theory and sim- 
ulation is quite good. The middle panel simply reflects the 
initial ui^-xx, i.e., it is the Fourier transform of equation (j68|l 
(with no jx = component, as described in Figure [51). 
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Fig. 10. 



Contour Plots at Two Times: Dotted lines 



trace contours of constant 



'-^swmg 



from a simula- 



tion identical to the one described in Figure [9j but with 
Ux = 600 instead of 150. To compare with the theoretical 
prediction (eq. | 65 | ■ or equivalently, eq. | 67 p. we plot two 
points on each contour where the contour's x-position is 
extremal. From equation (I67D . the extrema should be at 
^ext = {4> — tB,n~^ (x/a)) / Ky , where (j) = ""/S and 3n/2. We plot 
these two j/cxt as solid lines, with a = 0.082 (as obtained 
from Figure [SJ. The lines nearly go through the points, 
implying good agreement. 

6. SHEAR INSTABILITY IN REAL SPACE 

We seek a better understanding of the instability dis- 
cussed in ^ The instabihty of shear flows with inflec- 
tion points in the streamwise velocity field — or equiv- 
alently, with duJa,^\ I d x = (eq. [62j ) — is well known 
(jPrazin fc Reidl[20M ). However, we have not found in the 
literature a simple yet quantitatively accurate explanation 
for why it occurs. Why are only certain shear flows un- 
stable, and why only for sufflciently small wavenumbers? 
We give a dynamical explanation for a top-hat profile in 
%.![ and an explanation based on momentum and energy 
considerations in i j6.2l We also simulate the development 
of the instability into the highly nonlinear regime in i j6.3l 
We generalize these results from a top-hat profile to an 
arbitrary one in H6.41 and present another simulation in 



6.1. Top-Hat Vorticity Profile 

We consider in some detail the dynamics when the un- 
perturbed CO is given by the top-hat profile 



.'unp 







|a;| < s 
\x\ > s 



(75) 



where fj, and s are constants, with s > 0. The total unper- 
turbed vorticity is —q + tOunpi and we assume that <C q. 
The top-hat profile is perhaps the simplest profile that ex- 
hibits instability. It is not a realistic profile, both because 
real fluids are not discontinuous, and because it produces 
a My-field that is infinite in extent. But once the results 



for the top-hat are in hand, it trivial to extend them to 
profiles that do not suffer from these defects. We do so in 
^ ^6.4j|6.5l where we show that more realistic profiles be- 
have similarly to the top-hat. 
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Fig. 11. — Perturbed Top-Hat Profile: Main panel shows 
contours of constant vorticity. Dotted lines are for the un- 
perturbed profile (eq. [75]), and solid lines are perturbed. 
Top and right panels show the velocity fields that advect 
the contours. The top panel shows the amplitude of the 
background shear profile —qxy, and the right one shows Ux 
as a function of y. As long as inequality (I78D is satisfied, 
Ux can be taken to be independent of x. Uy is not depicted 
because advection in the y direction is dominated by the 
background shear. The data for this figure come from the 
simulation described in 96.31 at time=4000. This simula- 
tion has < and \ky\ < \fj.\/qs, so that instability results. 
Spatial and velocity axes are not to scale. 

We perturb the top-hat by displacing the two contours 
at X- and x+ by £,~{y,t) and £,'^{y-,t), as in Figure [TTl 
and seek equati ons of motio n for the contours ("contour 
dynamics," e.g., lPullinlll992D . We can adopt this partially 
Lagrangian description as long as a contour does not fold 
back upon itself; otherwise, ^ becomes multiple-valued, 
and a fully Lagrangian or an Eulerian description would be 
preferable. Although tracking contours is simpler than the 
mode coupling method of earlier sections, mode coupling 
can readily be generalized to three dimensions, whereas 
contour dynamics cannot: by tracking the contours, we 
rely on the fact that u is locally conserved, which is only 
true in two dimensions. We still assume that the dynam- 
ics occur in a box that is periodic in y with size Ly] but 
we now assume Lx is sufficiently large that the boundary 
conditions in x are unimportant. (From eq. [81| below, it 
suffices to take ^ l/fcy for all ky.) We shall make use 
of Fourier transforms, but only in the y-dimension, with 

uj{x,y)e-^^^y^ 



(76) 
(77) 



and similarly for Ux,ky{,x\ £,ky etc., where ky = (27r/Ly)x 
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integer. 

We shall see that instability can only occur when |fcj, |s < 
<gC 1; therefore we assume at the outset that pertur- 
bations are elongated in the streamwise direction, 

|fcy|s<l. (78) 

The contours are advected by the total velocity field 
~qxy + u (eq. [S]), implying 

[dt T qsdy)C^ = , (79) 

after dropping the nonlinear terms {—q^^ + Uy)dy£,^. In 
Fourier-space, 

d 
.dt 



signs, the flow is stable for all ky. When /i < 0, the flow 
is unstable provided that 



< K = 



qs 



(85) 

\fJ,\/2, which occurs 
or 



(80) 



We need not distinguish the Ux that advects f + fi'om the 
one that advects i^^, as we presently show. 



To find Ux^ky 



we have = ~c?yV '^uj (eq. [5]), or 



-sgr 



n(ky) / u;kAx')e-\''-^---'^\dx' (81) 



)sgn(fc;,) ,for |.t| < s (82) 



To derive the latter relation, the vorticity near the "step" 
at a; = s may be expanded as w = const — ^7?(x — s — ^"'") ~ 
fiS,^{y,t)S{x — s) + f{x) where H{x) is Heaviside's step 
function, 5{x) = dH/dx is Dirac's delta function, and f{x) 
denotes y-independent terms that may be discarded be- 
cause sgn(O) = 0. The step at a; = — s contributes with 
the opposite sign. The exponential in equation (|8ip may 
be dropped because of inequality ([TS]) . 

Even though uj varies rapidly with x, the Ux that it 
induces is independent of x. (More precisely, decays 
exponentially in x on lengthscale l/|Aj,| 3> s.) This result, 
which will simplify the analysis, corresponds to our earlier 
finding that swinging waves only contribute significantly 
when their phasefronts are radial. More generally, one ob- 
serves from equation (|8ip that any vorticity distribution 
that is elongated in the streamwise direction induces a ve- 
locity field that is much smoother in the x-direction than 
is Lo. The convolution smooths out variations in x with 

a smoothing length of l/\ky\; hence varies in a;— as words, if \ky\ > \fJ.\/qs, the flow is stable, as found above 



The maximum growth rate is 7,, 
at \ky\ — k/2, and 7 falls to zero as either ky 

We may now begin to appreciate how the shear insta- 
bility works in real space. Depending on the phases of the 
sine waves ^'^ (y) and ^~ (y) , the velocity field Ux that they 
induce can tend to amplify them. Figure [11] shows a grow- 
ing mode in an unstable shear flow that has /i < 0, so that 
the vorticity in the strip down the middle has the same 
sign as the background vorticity —q. The strip bulges in 
its central region, where |y| < Ly/A, and because of the 
sign of w, the bulge tends to produce a swirling velocity 
around itself with the same sense as the background shear. 
Therefore it produces a positive Ux in the top of the figure, 
and a negative one at the bottom, as depicted in the right 
panel. Furthermore, the strip is tilted into the top-right 
and bottom-left quadrants of the x — y plane. Although 
the tilt might not be obvious at first glance, if one com- 
pares the vorticity along any line of fixed \y\ with that at 
— |y|, one observes that the former is to the right of the 
latter. Because the induced Ux pulls the top half of the 
strip further to the right and the bottom half further to 
the left, it tends to amplify the tilt. In sum, Ux converts 
a bulge into a tilt. When < 0, the tilt is into the shear, 
having the same sense as a leading wave. 

To complete the picture, consider the effect of the back- 
ground shear, which carries ^+ downwards at speed gs, 
and ^~ upwards at the same speed. As long as the vor- 
ticity strip is tilted into the shear, the background shear 
converts a tilt into a bulge, thus completing the loop and 
allowing for exponential growth. 

Based on the above picture, we may estimate the fastest 
growth time as 7i^ix ~ £,^/ux ^ 1/ImI (eq- El]). Ad- 
vection of the two contours by the background shear can 
dephase them before growth can occur. Since the dephas- 
ing time is idcphasc ^/\Qsky\, it will be shorter than the 
fastest growth time if \ky\ is sufficiently large. In other 



well as in y — on the scale l/|fcy|. One can further argue 
on dimensional grounds that, in general, elongated vortic- 
ity distributions produce a velocity field with amplitude 
Ux ~ t-^S,- When ^ = 0, one must have m^, = by symme- 
try considerations. And although there is a dimensionless 
combination kyS (where s would represent the lengthscale 
of UJ in the x-direction), Ux cannot depend on s because it 
smooths out UJ on scales much larger than s. 

The simple equations ([50)1 and embody all the dy- 
namics of the initial stages of the shear instability. To 
solve them, we insert the trial solution cx e'''*, which 
produces for the eigenvector and eigenvalue 



1 



7 =F iqsky 
^— (qskyf - qs\ky\fj. 



(83) 

"f" ^ -[qslty)" - qs\ky\fj. (84) 

The behavior of the growth rate 7 is very similar to that 
found for previously considered vorticity profiles (eqs. [JS] , 
[73] and Fig. [8]). When n > 0, i.e. when the perturbed 
vorticity and the background vorticity —q have opposite 



6.2. Momentum and Energy for the Top-Hat 

It is instructive to see how momentum and energy inter- 
act to trigger the shear instability. Not only does this give 
a deeper understanding of why the instability operates, it 
also gives a sense of the nonlinear outcome of the insta- 
bility, since it quantifies the driving force behind it. Fur- 
thermore, the interaction between momentum and energy 
is generally important in accretion disks. Understanding 
how this works in a nontrivial incompressible flow might 
be helpful when considering more complicated effects, such 
as three dimensional motions, baroclinicity, or magnetohy- 
drodynamics. 

We first derive expressions for the momentum and en- 
ergy in terms of the dynamical variables , and show that 
the equations of motion ([50)1 and ([5^ exactly conserve 
these quantities. Perturbations to the top-hat profile have 
y-momentum (per unit Ly and setting p = 1) 



M— {uy — Uy\^±^Q)dxdy / L 



(86) 
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{Uy 



=o)dx 



(87) 



where the bar denotes a y-average, 

u{x) = / u{x,y)dy/Ly 



Since Ux = (because Ux — —dyV^'^uj), the total x- 
momentum always vanishes, and need not be considered 
further. Henceforth we call M the momentum. From 
Co — dxUy, we may express 



Uy-Uy\^±^i^^ I {uj{x') - uj{x')\^±^Q)dx' (89) 

^i^{\£,uf^{^ + ^)-K?s{x-s))m 

where the latter relation applies for a single ky component 
of (including the —ky piece). To derive it, note that 
the integrand in equation (|89p vanishes everywhere except 
where a line of constant x' pierces either of the two con- 
tours. Considering first the vicinity of the ^+ contour, we 
have ZJ(a;') — IZJ(x')|j±^o = ~tA^k ^^^5 / dx\x' -s, after writ- 
ing uj in terms of Heaviside's function as below equation 
(j82p . and then expanding. The contour contributes the 
analogous amount, but with the opposite sign because the 
jump in vorticity across it is in the opposite sense. There- 
fore both steps contribute 

M = M_ + M+ (91) 

M±^TA^I4±r (92) 

The overall signs of Af_ and Af+ will play an important 
role in what follows. They are determined solely by the 
sign of the jump in vorticity across the corresponding step, 
independent of the functions ^~ and We may under- 
stand this is as follows. To be definite, we consider first 
the step at a; = — s and set /i > 0. Because a) = duy/dx, 
the unperturbed Uy is constant for x < —s and has posi- 
tive slope (= ^) for X > —s. Any perturbation must 
widen the range in x over which u) makes a transition from 
to pL. Hence it will widen the range over which Uy in- 
creases in slope. Therefore the perturbed Uy must exceed 
its unperturbed value and M_ > 0. Generalizing this rea- 
soning, it is clear that wherever the unperturbed profile 
has a positive jump in vorticity, perturbations always con- 
tribute positive momentum; similarly, perturbed negative 
jumps contribute negative momentum. 
The perturbed kinetic energy is 

j i-Qxy + uf - {-qxy + uf\^±=Qdxdy (93) 

= E-+E++ (94) 



where 



E,= 



'qx{uy 



--=o)dx 



which follows from equation (|90p . and 

E.,,2 = / ^(u2 - u'^\^±^Q)dx 

,2 



2,1 kq, 



(95) 
(96) 



(97) 



(98) 



which follows from Uy^ky = {i/ky)dux.ky/dx and equation 
(|82|) (including now the exponential previously dropped 
from equation ((8T|) ) . There are two different kinds of en- 
ergy. The first, composed of E_ and is localized near 
each step in w, and is due to changes in Uy that are second 
order in we call it the mean-flow energy. The second, 
> 0, is spread out in x over a range '-^ ^/\ky \ ^ s, and 
is due to the terms in u that are first order in we call 
it the perturbation energy. 

The equations for conserve momentum, 

dM± , , 

* ±F (99) 



dt 



where 



F = sgn(fcy)-^' 



(100) 



is the rate at which momentum flows from the step at 
x — —s to the one at +s. In terms of more familiar vari- 
ables, it can be shown that the x-dependent momentum 



flux 



F 




\x\ < s 
\x\ > s 



(101) 



which is in accord with equation (j99p . Similarly, the equa- 
tions for conserve energy, 
d 



dt 



E~ 



-q{2s)F 
dE„2 



dt 



(102) 
(103) 



which may be understood as follows. Consider two stream- 
lines of the unperturbed flow, one a,t x = xi and the other 
at X = a;2, with xi < X2, and label their unperturbed 
y- velocities Vi and V2- Their unperturbed specific ener- 
gies (energy per unit mass) are then v'f/2 and u|/2. We 
perturb the streamlines by transferring y-momcntum from 
one streamline to the other. If the specific momentum 
transferred from xi to X2 is Aw, then the change in mean- 
flow energy is 

(wi - Avf - vj {v2 + Av)^ - 



° Note that outward transport of angular momentum in an accretion 
lower the relative velocity of two streamlines. 



2 2 

Ki{v2-vi)Av^-q{x2-xi)Av , (104) 

which is equivalent to equation (|102p . When the momen- 
tum transfer tends to lower the relative velocity of the two 
streamlines (i.e., Av > 0), energy is extracted from the 
mean flow. Conversely, when the transfer tends to increase 
the relative velocity, energy must be added to the mean 
flowH Of course, the total energy is conserved: if energy 
is extracted from the mean flow it must correspondingly 
increase E^i. 

Let us see what the above results imply for the shear 
instability, first setting /i < so that the instability can 
operate. Perturbations to the step at x = —s have negative 
momentum, and those at x = s have positive momentum 
(eq. [22] )• Therefore perturbations transfer momentum 
from X = —s to X = s. Since this transfer tends to lower 
the relative velocity of the streamlines at x = ±s, en- 
ergy is extracted from the mean flow and added to £'„2. 
This is the recipe for instability. However, for instabil- 
ity to actually occur the increase in i?K2 must balance the 

disk is equivalent to a momentum transfer in a shear flow that tends to 
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decrease in mean-flow energy. Since Ey^2 cx l/|fcj,|, this 
is only possible for sufRciently small \ky\. Quantitatively, 
unstable perturbations must have M_ + M_|_ = 0; other- 
wise, M_ -I- M4. could not remain constant as | grows. 

Therefore, j^^T | = \- Similarly, energy conservation 
implies i?„2 + + = 0, or 



2qs 



(1 



(105) 



this 



and 4^ 
. The largest 
|/z|/gs, as for 



where (p is the phase difference between 
relation also follows from equations (|83l) - ((84l) 
\ky \ occurs when (f> = n, and is equal to k = 
equation ([55]) . 

If /i > 0, equation (|105p can never be satisfied, and 
there is always stabil ity, as implied by Fj0rtoft's theorem 
(jPrazin &: ReidI [20041 ) . This seemingly mysterious behav- 
ior can also be understood from momentum and energy 
considerations. When /i > 0, the perturbed momentum at 
X — —s is positive, and that at x = s is negative. Therefore 
any perturbation must transfer momentum from x = s to 
X — ~s, which enhances the background shear, and must 
increase the mean flow energy E- + E^ . By energy con- 
servation this is only possible if the perturbation energy 
i5„2 decreases, and since £'„2 is always positive, no pertur- 
bation can grow. 

6.3. Numerical Simulation of the Top- Hat 

Figure[T3shows the nonlinear evolution when the unper- 
turbed vorticity is given by equation (|75p . The initial per- 
turbation to the curves at a; = ±s is unstable {ky = 0.75k), 
and wraps up into a vortex. 

To make this flgure we wrote a Lagrangian code, which 
can handle the rapid variations in oj more easily than can 
the pseudospectral code. The Lagrangian form of the 
equation of motion (eq. [5]) is dtx = —qxy + u, where 
X = x{xo,t) is now a function of the initial coordinate 
and time, dt is taken at fixed Xq, and u is determined by 
equation ([9]). We chose 2000 initial points along each of 
the two curves at x = ±s, and evolved these points assum- 
ing periodic boundary conditions in y and open boundary 
conditions in x. Given x for each of these points, we im- 
mediately know w, since uj is locally conserved. To convert 
uj to u, we evaluated equation ^ by interpolating uj onto 
a grid, taking the Fourier transform of uj, multiplying by 
the appropriate factors of fc, taking the inverse Fourier 
transform to yield u on the grid, and finally interpolating 
from the grid back onto the curves. To check the code, we 
plot in Figure [12] at time= 3000 the eigenfunction 



^^(y) cx ±cos \ky\y^tan 



1 



qs\ky\ 



(106) 



(eq. [83]). These curves are indistinguishable from the 
code's output, implying agreement between code and the- 
ory. We have also checked that the growth rate at early 
times, while < s, is oc exp(7t), where 7 = 1/1200 
(eq. [84j). As usual, for our numerical results we measure 
time in units of , i.e., we set q = 3/2. 

At late times, uj exhibits rapid spatial variations, and 
the boundaries that separate the regions where w = 
from the ones where uj — fi become more and more convo- 
luted. The Lagrangian code has no viscosity; the evolution 
shown in Figure [T^ breaks down at time > 10, 000 when 



the curves are so convoluted that adjacent points on a 
curve are widely separated. If we were to include an ex- 
plicit viscosity, it would wipe out the rapid variations in 
UJ at late times, leaving a smooth vortex. At what stage 
the viscosity acts depends on how small the viscosity is. 
In astrophysical disks, the viscosity is extremely small, so 
if a shear instability acts, it might lead to extremely rapid 
variations in uj. Whether this is unstable to three dimen- 
sional perturbations is an interesting possibility. 



time = 3000 



time=6000 time=8000 



time = 9000 



-0.002 






Fig. 12. — Numerical Simulation of the Top- Hat: The un- 
perturbed vorticity is given by equation (I75D . with s = 0.1 
and fi = —0.002; the initial perturbation to the curves at 
X = is, ^^{y), are sinusoidal with ky = 0.01 = 0.75k, and am- 
plitudes = 0.001. The initial perturbation is unstable to the 
shear instability, and wraps up into a vortex. The numeri- 
cal simulation was done with a Lagrangian code described 
in the main text. 

6.4. Infinite Number of Steps 

Consider now the general case, in which the unperturbed 
UJ is given by Wunp = , for Xi < x < x^+i , where 
and Xi are sets of constants with i = 1, • • • , 00. As before, 
there is a curve Ciu^ t) at each step Xi, whose equation of 
motion is, generalizing equation (|80p . 

^ - qx.ik^ Cky = u.^,ky ■ (107) 
Equation (|82p generalizes to 



(108) 

(109) 



where the jump in vorticity at step i is 

ujI=uj'- cj'-i 
We continue to assume that |A:j,|s <C 1, where s is now the 
extent of the region of non-zero vorticity. To find insta- 
bility, we set Q. oc e'*'*, yielding for the eigenfunction and 
eigenvalue 

1 



7 - qxitlty 
-^agn{ky)Y^ 



7 - qxiiky 



(110) 
(111) 
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thus recovering the results derived with swinging waves 
(eqs. [67], [64]) when we make the identification 7 = aqKy. 

6.5. Top-Hat With Wings 

Because the top-hat vorticity profile has / ojdxdy 7^ 0, 
it produces a velocity field Uy that is infinite in extent. Al- 
though this is somewhat disconcerting, profiles that do not 
suffer from this problem behave similarly to the top- hat. 
To illustrate, we consider the top- hat with wings, i.e., 

[ A* , < s 
Wunp = S ' s <\x\<2s , (112) 

{ , 2s < 
which clearly integrates to zero. Figure [T3l shows a numer- 
ical simulation of this profile, done with the Lagrangian 
code described in i j6.3l See caption for details. For unsta- 
ble initial conditions, the vorticity again tends to wrap up 
into a vortex. 

It is perhaps not surprising that the infinite extent of the 
unperturbed Uy is not important, since the unperturbed 
Uy was not important for any of our three explanations of 
the shear instability, in ijifSl 16.11 and 16.21 



time = 3000 time = 4000 tinie = 6000 




Fig. 13. — Numerical Simulation of Top-Hat with Wings: 
The unperturbed vorticity is given by equation (I112II ■ with 
s = 0.1 and /i = —0.002; the initial perturbations are sinu- 
soidal with ky = 0.01. The dispersion relation (eq. 
shows that this situation is unstable, with growth rate 
7 = 1/810. We checked that at early times the numeri- 
cal simulation reproduces this growth rate. As in Figure 
1121 the nonlinear outcome is a tightly wrapped vortex with 
convoluted boundaries. 

7. CONCLUSION 

Two dimensional shear flows are linearly stable, but 
nonlinearly unstable. Linear modes that are aligned with 
the shear ( "axisymmetric" ) have zero frequency; those 
that are not aligned ( "swinging" ) are transiently amplified 
before decaying. When the two couple nonlinearly, they 
are always unstable provided that the swinging mode's az- 
imuthal wavelength is long enough. If this is the case, the 



swinging mode couples with the axisymmetric one to pro- 
duce a new leading swinging mode with larger amplitude 
than itself, implying instability. Quantitatively, the crite- 
rion for instability is 

fcy,sw < M ^ 1^1 

kx.axi 1 3f2/2 

where q is the ratio of perturbed to background vortic- 
ity, and fcy,sw/fcx,axi is the ratio of swinging-azimuthal to 
radial-axisymmetric wavenumbers. This shows that even 
when the amplitudes of the perturbations are vanishingly 
small, instability results for small enough ky^sw 

When the instability is triggered, the swinging waves 
grow exponentially at the rate ~ , assuming that their 
wavelength is larger than the critical value (eq. |113p by 
a factor of order unity. Eventually, regions of negative 
vorticity concentrate into vortices with increasingly con- 
voluted boundaries. Since astrophysical disks have very 
little viscosity, such vortices might become turbulent in 
three dimensions before viscosity can act to smooth them. 

In fJU we explained how the instability works in real 
space. When a nearly axisymmetric strip of perturbed vor- 
ticity has a small bulge, this bulge tends to tilt the strip. 
Advection of the tilted strip by the background shear re- 
inforces the bulge when equation ()113|1 holds, where 
is now interpreted as the width of the strip, and ky^sw is 
the azimuthal wavenumber of a perturbation to the strip. 
(It is also required that w < if the strip only has a sin- 
gle sign of vorticity.) We also showed that the instability 
can be understood from momentum and energy considera- 
tions: when equation (jll3p holds and the strip is unstable, 
the momentum transferred from a streamline on one side 
of the strip to the other tends to lower the relative veloc- 
ity of two streamlines (which is equivalent to transferring 
angular momentum outwards in an accretion disk), and 
therefore converts shear energy into perturbation energy. 

It is tempting to speculate that this nonlinear instabil- 
ity might provide a route to turbulence. Perhaps vortices 
are unstable to three dimensional perturbations, and these 
perturbations can create new axisymmetric modes that are 
in turn unstable to the production of vortices. It is also 
possible that other effects are required. Perhaps vorticity 
is produced by convection or by stirring from planets, and 
if enough vorticity is produced, the disk becomes unstable 
in the manner described in this paper. Even if other ef- 
fects are required, we feel that the technique introduced in 
this paper — the examination of couplings between swing- 
ing modes — will be helpful. It gives a well-defined pro- 
cedure for analyzing nonlinear effects, both theoretically 
and numerically. Similar techniques might also be em- 
ployed to give insight into other issues in shear flows, such 
as the nonlinear development of the magnetorotational in- 
stability. 

We thank Gordon Ogilvie for showing us how easy it is 
to implement shearing box boundary conditions in a pseu- 
dospectral code, and Jeremy Goodman for helpful discus- 
sions at an early stage of this project. We also thank the 
referee for helpful suggestions. 



15 



APPENDIX 

DESCRIPTION OF THE PSEUDOSPECTRAL CODE 

The pseudospectral code simulates the 3-D shearing box equations of motion (eqs. [5]-[6]), though in the present paper 
we only present 2-D simulations. Most of our methods are standard, except for our new method for remapping modes 
that is both simpler and more accurate than that used by other investigators. 

In our simulations, the unit of time is chosen so that 

= 1 , (Al) 

and we also set the constant density 

p=l. (A2) 
Equations (O-® are solved in a box of size ^ Ly x that is subject to periodic boundary conditions in y and z, 

u{x,0, z) = u{x, Ly, z) (A-3) 
u{x,y,0) = u{x,y,L^) (A4) 
and to a shifted periodic ("shearing sheet") boundary condition in x: 

u{0,y,z) ^ u{L.j;,{y - qtL^)modLy,z) . (A5) 
The modulus function "mod" is defined as follows: 

{y - qtL^) mod Ly = y - qtL^ + pLy (A6) 
where p is the integer that makes the right-hand side of the above expression lie within the box in the j/-dimension: 

<y - qtLx +pLy < Ly. 

Changing to Lagrangian coordinates X,T that shear with the background velocity —qxy (eqs. [I2]-[I3])j equation ([5]) 
becomes 

du 

— = 2uyX - (2 - q)uxy - VxP - u ■ V^u , (A7) 

where Va, = V is a gradient with respect to x, not X. 

Lagrangian coordinates {X, Z) lie within a box of size x Ly x L^] the boundary conditions on the sides of this 
box are periodic in all three dimensions. The Fourier transform of u with respect to X is 

uj^,Jy.J. = r T T I «cxp(-zK • X)£X , (A8) 

where the wavevector K takes on the discrete values 

27r 27r 27r , 

= — , Ky = —Jy , = —J, , (A9) 

Lx J-iy L^z 

with JxT-Jyi Jz integers. Hence the Fourier transform of equation (jA7|) with respect to X is 

du - ' 

— = 2uyX - (2 - q)uxy ~ ik{K, T)P - ik{K, T)-[uu] , (AlO) 

where, for clarity, we denote Fourier transforms with a tilde instead of the subscripts Jx, Jy, Jz- The (x, y, z) components 
of the wavevector k{K , T) are defined as follows 

k{K, T) = {Kx + qTKy, Ky, Kz) . (All) 

The last term in equation (jAlOp has been written in a tensor form; it represents a vector whose j'th component is 

3 

- i'^kiuluj , (A12) 
1=1 

where uiUj is the Fourier transform of the product uiUj. Equation (|A10p is supplemented by the Fourier-space version of 
equation ([S]): 

k{K,T)-u = {) (A13) 

The numerical code evolves equation (|A10|) in time using second-order Runge-Kutta, with P determined by the incom- 
pressibility constraint (eq. |A13j V Equation (jAlOp is integrated for UxXUy x Uz different modes, each of which is labelled 
by the integers {Jx, Jy, Jz), with —nx/2 < Jx < nx/2, —ny /2 < Jy < ny/2, and —nz/2 < Jz < nz/2^ 

The code equations differ from the exact equations (jA8[) - (jA13|) in two ways. First, the Fourier transform in the code is 
the discretized version of equation (|A8[) : 



'f^x^y^z 

^ In truth, only half of these modes are simulated, which makes the code run twice as fast as it otherwise would. The other half that are not 
simulated are determined by the condition that the components of u are real, not complex, numbers. But for the purposes of the discussion in 
the present section, we ignore this complication. 
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with the sum taken over Ux x riy x Uz discrete Lagrangian points X. And second, kx is not given + qTKy, as would be 
inferred from equation (lAllD. I n this respect, our code differs from other codes (e.g., lRogallol[l98lt lUmurhan fc RegevI 
arranco fc Marcudl200l . A single mode with a given {J^, Jy, Jz) has the following spatial dependence 



exp 



27ri 



JxX 



JyY 



(A15) 



exp 



27ri 



Jy 



-y 



Lz 



(AI6) 



Jx 

^Lx ±Jy ^ 

In substituting for Y with equation (|12p. the modulus term may be dropped because it does not contribute to the 
exponential. By the same token, we can replace Jx on the right-hand side of the above expression with Jx +pnx, where 
p is any integer, without affecting the exponential. (Note that x — {0,1, ...,nx — l}Lx/nx.) This is the well-known 
phenomenon of aliasing. So instead of equation (|A11|) for kx, we can equally well use a similar expression, but with 
Jx Jx + pux- The value of p should be chosen so that the physical wavenumber kx is shifted as close as possible to 
zero. Other p's correspond to lengthscales smaller than the grid's scale. Hence, 



kx — 



27r 

Lx 



Jx + qTJi 



Lx 



mod 



(A17) 



where mod is a modulus operator: 



J mod nx = J + pUx 



(A18) 

with p an integer such that —Ux/'i < J + pUx < nx/2. (These limits differ from those used to define the earlier modulus 
operator below eq. |A6| .) Since other authors use kx — Kx + qTKy instead of equation (|A17|) for kx, the kx of each of 
their modes grows to large positive or negative values. To stop this growth, they periodically remap the modes to small 
kx ■ But this has the undesirable consequence that the distribution of kx 's that is being used to simulate the fluid changes 
substantially from just after a remapping event to just before the subsequent remapping. Conversely, with equation (|A17|) 
there is no need for an explicit remapping; in effect, the remapping is done automatically by the modulus operator. 

The nonlinear t erm in equation (lAlOD . written explicitly in equation (jA12[) . is evaluated using the standard pseudospec- 
tral method (e.g., iMaron fc Goldreich.2001i . and references therein): after dealiasing (see below), uj^^j^^j^ is transformed 
into u by performing an inverse fast Fourier transform (FFT). Then, the product UiUj is formed in real space, and finally 
the FFT of this product is taken. 

To prevent aliasing errors when evaluating the nonlinear term, mj^.j,,,,/^ is truncated just before the time derivative is 
evaluated by setting to zero any mode whose \kx\, \ky\, or \kz \ exceeds 2/3 of its maximum value — the standard "2/3-rule"; 



explicitly, it 



wherever one of the following hold 



\Jy\>f 
\Jz\ > ^ 



Jx + qTJy 



mod 



> 



(A19) 
(A20) 

(A21) 



The last condition above follows from equation (jA17p . Note that the kx of a mode is remapped by the modulus operator 

from +kx to —kx when {Jx + qTJyLx/ Ly) mod Ux — rix/^. Such modes have had their amplitude set to zero because 

they satisfy the inequality of equation (jA21[) . Hence there is no danger that remapping artificially introduces power into 
leading modes. 
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